AN  AGGREGATION-BASED  DOMAIN  DECOMPOSITION  PRECONDITIONER  FOR 

GROUNDWATER  FLOW  * 
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Abstract.  In  this  paper  we  give  a  convergence  analysis  of  a  two-level  additive  Schwarz  method  in  which  the  coarse 
mesh  basis  is  constructed  with  aggregation,  a  method  common  in  algebraic  multigrid.  This  coarse  mesh  does  not  need 
geometric  information  about  the  subdomains  and  can  readily  be  used  on  unstructured  spatial  meshes.  We  illustrate  the 
performance  with  several  computational  examples. 
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1.  Introduction.  This  paper  provides  a  convergence  theory  for  the  two-level  Schwarz  pre¬ 
conditioner  first  described  in  [8-10].  The  preconditioner  is  a  two-level  additive  Schwarz  method 
[4, 5, 18].  Its  novel  feature  is  the  use  of  aggregation  ideas  from  algebraic  multigrid  [1, 22]  and  bal¬ 
ancing  Neumann-Neumann  methods  [13-15]  to  construct  the  coarse  mesh.  The  use  of  aggregation 
arose  from  necessity.  In  the  applications  reported  in  [8-10]  the  subdomains  were  irregular,  and  a 
coarse  mesh  based  on  “hat  functions”  over  the  subdomains  was  impractical.  For  the  same  reason, 
we  needed  minimal  overlap  between  subdomains.  Unlike  the  method  from  [3],  we  do  not  need  to 
create  a  coarse  mesh  geometry  or  use  geometric  information  about  the  subdomains. 

While  our  theoretical  results  allow  for  overlap,  and  some  of  the  experimental  results  in  §  3  are 
for  subdomains  with  substantial  overlap,  the  application  of  the  preconditioner  in  practice  has  been 
for  minimally  overlapping  subdomains. 

Our  analysis  uses  the  standard  finite  element  framework  from  [18,24],  The  preconditioner  also 
works  well  in  the  context  of  finite  differences,  however,  as  some  of  the  examples  in  §  3  illustrate. 

This  preconditioner  was  developed  for  use  in  the  Adaptive  Hydrology  model  (ADH)  [19],  a 
production  code  being  developed  at  the  US  Army  Engineer  Research  and  Development  Center. 
Because  ADH  is  a  three-dimensional  unstructured  mesh  code,  detailed  grid  refinement  studies  are 
too  costly  to  perform;  therefore,  in  this  paper  we  illustrate  the  performance  of  the  preconditioner 
with  computational  examples  in  two  space  dimensions,  a  setting  in  which  grid-refinement  studies 
and  investigations  of  the  effects  of  overlap  can  more  readily  be  done. 

1.1.  Problem  Formulation.  We  begin  with  the  weak  form  of  an  elliptic  boundary  value  prob¬ 
lem  on  a  domain  Sic  Rd  with  boundary  T.  The  goal  is  to  find  u  E  V  such  that 

(1.1)  a(u,  v)  =  l(v)  for  all  v  G  V 
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where  l  is  a  linear  functional  on  V,  and  V  is  an  appropriate  function  space.  In  our  examples  the 
problems  will  be  second  order  with  Dirichlet  or  Neumann  boundary  conditions  and  V  C  H1( Q). 
We  take  approximating  spaces  of  continuous  functions  Vh  C  V.  The  approximating  problem  at 
level  h  is  to  find  uh  E  Vh  such  that 

(1.2)  a(uh,  v )  =  l(v )  for  all  v  E  Vh. 

The  problem,  (1.2),  is  equivalent  to  a  linear  system 

(1.3)  Auh  =  / 

on  Vh  where  a(u,  v )  =  (Au,  v )  for  all  u,  v  E  Vh.  Here  (-,  •)  is  the  L2(Q)  inner  product. 

Schwarz  preconditioners  are  designed  to  accelerate  Krylov  space  iterative  methods  for  the 
solution  of  (1.3). 

1.2.  One-level  Method.  We  begin  with  the  one-level  additive  Schwarz  preconditioner.  We 

divide  Q  into  subdomains  with  an  overlap  width  of  at  least  5  and  assume  that  U/=i  %  =  H. 

Let  Rj  be  the  restriction  map  from  an  element  of  Vh  to  the  subspace  Vj  of  Vh  with  support  on 
Q,j.  Let 

Aj  —  RjARj 

be  the  subdomain  operator.  We  assume  that  Aj  is  nonsingular  on  V3  and  define 

Bj  =  Rj  Aj  Rj , 

where  Aj  is  an  approximation  of  Aj.  The  one-level  additive  Schwarz  preconditioner  is 

(1.4)  M  = 

3=  1 

1.3.  Two-level  Method.  The  two-level  additive  Schwarz  method  adds  a  coarse  mesh  term 

Bq  =  RqA0  1i?0 

to  the  one-level  preconditioner.  Here  H0  is  an  approximation  of  A0.  We  let  VH  denote  the  space  of 
coarse  mesh  basis  functions.  If  the  coarse  mesh  restriction  map  i?0  and  the  coarse  mesh  operator 
A0  are  well  designed,  the  condition  number  of  MA  is  significantly  reduced. 

One  way  to  define  a  coarse  mesh  problem  is  to  discretize  the  continuous  problem  on  a  coarser 
mesh.  There  are  a  few  difficulties  associated  with  forming  the  coarse  problem  this  way.  For 
unstructured  meshes,  such  as  the  ones  considered  in  [8-10],  the  interpolation  operators  between 
the  fine  mesh  and  the  coarse  mesh  are  difficult  to  define.  Second,  a  coarse  mesh  must  be  generated, 
stored,  and  parallelized.  Finally,  the  PDE  must  be  discretized  and  recomputed  on  the  coarse  mesh. 

Alternatively,  the  discretization  of  the  coarse  mesh  operator  may  be  defined  in  terms  of  the 
existing  fine  mesh  discretization.  A  Galerkin  or  variational  coarse  grid  correction  uses  the  fine  grid 
matrix  to  obtain  the  coarse  grid  operator  as  A0  =  R0AR where  R0  is  the  interpolation  operator  to 
the  coarse  mesh  function  space,  and  R q  is  the  restriction  operator.  Methods  that  obtain  the  coarse 
mesh  matrix  by  using  the  fine  mesh  matrix  are  called  nested  or  multilevel  methods  [18].  If  the 
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coarse  mesh  basis  functions  are  obtained  from  the  fine  mesh  basis  functions,  then  the  coarse  mesh 
space  VH  is  contained  in  the  fine  grid  space  Vh. 

We  use  the  aggregation-based  basis  from  [8-10]  in  this  paper,  where  one  coarse  mesh  ba¬ 
sis  function  is  defined  for  each  subdomain  as  the  sum  of  the  fine  mesh  basis  functions  for  that 
subdomain. 

To  set  the  notation  that  we  will  need  in  §  2,  let  the  expansion  of  a  function  u  G  Vh  in  the  finite 
element  basis  be 

(1.5)  u  =  YJuii)i 

i 


where  the  tf's  are  the  basis  functions  for  the  fine  mesh,  and  a  function  uc  G  VH  can  be  represented 
on  the  coarse  mesh  space  as 

(1-6)  uc  =  Y,ucK^k 

K 

where  the  \b/r’s  are  the  finite  element  basis  functions  for  the  coarse  mesh  space.  Since  V 11  C  Vh, 
can  be  written  as 

(1.7)  v1;k  =  y"  Rk^i. 

j 

The  index  K  represents  the  subdomain  number.  For  our  coarse  mesh  basis  functions,  the  value  of 
Rkj  is  either  0  or  1.  If  the  fine  mesh  basis  function  xpj  is  contained  in  subdomain  K,  then  Rk3  =  1; 
otherwise,  Rxj  =  0. 

Further  expanding  the  representation  of  uc  gives 

Uc  =  J2kuCk^  k 

—  Sir  Uqk  Sy  RKj'ipj 
—  Sj  (Sir  ucKRi(j )  tpj 
=  S j  (RTUc),  Ipj 

Thus  RT  is  the  operator  which  interpolates  from  the  coarse  mesh  to  the  fine  mesh.  Any  function 
on  the  coarse  mesh  can  be  represented  solely  in  terms  of  the  already  existing  fine  mesh  functions, 
making  the  formulation  of  a  separate  coarse  mesh  unnecessary. 

1.4.  Condition  Estimate.  In  Assumption  1.1  we  make  precise  the  idea  that  H  is  the  charac¬ 
teristic  diameter  of  a  subdomain.  In  Assumption  1.2  we  make  precise  the  overlap  5  between  the 
subdomains  and  the  properties  of  the  coarse  mesh  basis  functions. 

Assumption  1.1. 

1.  There  is  C  >  0  such  that  diam(Qj)  <  CH  for  all  j  =  1. ... .  J. 

2.  There  is  c  >  0  such  that  for  all  x  <G  Q  there  exists  j  >  1  such  that  x  G  11,  and 


dist  ( x ,  dQj\dQ)  >  c5. 
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3.  There  are  Cr,  C\ .  C'2  >  0  such  that  for  all  x  G  Q  the  ball 


B  (x,  CrH )  =  {y  G  fl  :  dist  (y,  x)  <  CRH} 


intersects  at  most  C\  +  Cj  subdomains  Qj  ( i.  e. ,  an  object  of  diameter  O  (if)  intersects  at 
most  0(1)  subdomains  Qt). 

4.  p(Qj)>CHd,j  =  l,...,J. 


In  Assumption  1.1,  p  denotes  Lebesgue  measure. 

ASSUMPTION  1.2.  Assume  the  basis  functions  'It,  of  the  coarse  space  satisfy 


1. 


4/ 

II*. 


<  C 


Hd~ 


OiH(n)  —  ^  S 

2L2  <  CHd 


2.  There  is  a  domain  Qrnt  C  Q  such  that  Jj.,  'It,  (x)  =  1  for  every  x  G  Qvnl  and  dist (x .  dQ)  < 
C6  for  every  x  G  Q\Qint. 

3.  supp  (\&j)  C  fli. 

Parts  1  and  2  of  Assumption  1.2  differ  from  similar  assumptions  in  [1]  in  the  if1  semi-norm 
of  the  coarse  mesh  basis  functions,  and  in  the  distance  from  the  space  Qrnt  to  9Q. 

In  §  2  we  prove 


THEOREM  1.1.  Let  Vh  —  Vj  C  C(Q)  and  let  Assumptions  1.1  and  1.2  hold.  Assume 
that  there  is  u  >  1  such  that 


(v,  v )  <  (Aj  1Ajv,  v )  <  co(v,  v ) 
for  all  v  G  Vh  and  0  <  j  <  J.  Let 

(1.8)  M  =  J2Bj 

j= 0 

then  there  is  C  >  0,  independent  of  H  and  h  such  that 


(1.9) 


k(MA)  <  Cu(l  +  (H/h)2). 


2.  Convergence  Theory.  We  base  our  analysis  on  the  result  from  [23,24]. 

THEOREM  2.1.  Let  K0  be  a  positive  constant  so  that,  for  any  v  G  Vh,  there  exists  a  decom¬ 
position  v  =  E/_o  vi  such  that  Vi  G  V)  and 


(2.1)  (Aividi)  <  Ko  (Av,v) . 

i= 0 

Let 


(2.2) 


j 

Ki  =  max  V  c , , , 

djdt,  3 


where,  for  1  <  i,  j  <  J,  £ij  =  0  if  PiPj  —  0  ( i. e. ,  ifXjlAj),  Eij  =  1  otherwise.  Then 
(2.3)  k  ( MA )  <  coK0  (1  +  K1), 
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where 

u  =  max  Amax  (A^Aj)  . 

0 <j<J  V  J 

We  assume  that  the  energy  norm  is  equivalent  to  the  H 1  seminorm,  and  we  can  therefore  replace 
(2.1)  by 

j 

(2-4)  X  Wjjqn)  —  -^oM^qn)- 

«=o 

Our  estimate  for  K0  will  be  based  on  (2.4). 

The  value  of  K\  is  an  indicator  of  the  number  of  subdomains  that  contain  any  given  point  in  Q; 
we  assume  our  partition  is  such  that  K1  =  0(1).  We  solve  the  subdomain  problems  exactly  in  our 
numerical  results,  so  u  =  1  in  §  3.  Thus  our  condition  number  estimate  is  based  on  the  estimate 
of  Kq,  which  we  obtain  using  Lemmas  2.2  and  2.3.  Our  analysis  uses  some  of  the  techniques 
developed  in  [1],  but  we  do  not  smooth  the  coarse  mesh  basis.  Our  estimates  in  Lemma  2.2  are 
thus  different  from  those  in  [1], 

We  define  the  coarse  mesh  operator  Q  :  Vh  — >  VH  by 


f  1  f 

Qu  =  V  cq\Eq,  cq  =  a*  (u)  =  j-  .  /  u  ( x )  dx, 
~~  u  ( o,)  Jih 


i—  i  A*  (^t) 

where  u  G  Vh.  Using  the  Cauchy-Schwarz  inequality  and  Assumption  1.1  we  get 


/  u  ( x )  dx 


Hence 


ffd/2 

I  Oil  <  ll“l 


=  Hu I 


II2 


L2  ' 


The  value  of  K0  depends  on  the  bound  of  the  energy  of  Qu  and  on  the  L2  bound  of  the  error  in  the 
coarse  mesh  operator,  i.e.,  u  —  Qu.  These  bounds  are  provided  in  Lemma  2.2. 

Here,  as  in  the  remainder  of  this  section,  C  is  a  constant  that  is  independent  of  H  and  h.  C 
may  increase  as  the  analysis  progresses. 

LEMMA  2.2.  If  Assumptions  1.1  and  1.2  hold,  then 


(2.5) 

(2.6) 


| u  —  Qu \\2L2  <  CH2  \u\2h1 
\Qu\2Hi  <  Cy  \u\2H!  . 


Proof.  Let 


b  =  n\Qint 
B  =  {i  :  Qi  0  B  f  0} 

B'  =  u 

ieB 

W  =  sup  {dist  (. x ,  dfl)} 

x£B' 

B0  =  {x  G  O  :  dist  (x,  dQ)  <  W} 
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If  x  e  B,  then  (list  (x,  dQ)  <  C6  by  Assumption  1.2.  Thus 

B'  =  Uj  {0,  :  dtti  n  dtt  ±  0}  . 

Therefore,  we  have  that  distxeB'  (x,dQ)  <  CH  since  diam  ({1,1)  <  CH  for  all  i.  Hence  IT  <  CH 
and  using  the  Poincare  inequality  [6]  we  get 

(2-7)  <  IMIl2(b0)  —  CH  \u\fj((B0)  ■ 

Let  A fi  =  {]  :  flj  fl  Q,  0  }  and  let  card(Afi )  denote  the  cardinality  of  A'"; .  We  have 

WQuWIhb)  =  (Siee  on  (u)  (u) ,  Zjes  olj  (u)  ^  («)) 

—  E*eB  I ai  (M)l  \aj  (M)l  II^jIIl2 

<  §  E iets  Ej&winB  (al  (u)  +  aj  («))  CHd  since  a2  +  b2  >  2 ab 

<  CHd  max  {card(Ni)}  Ejgb  ai  (u)  <  C  Eieis  ll«lli2(n.)  <  C \\u\\2l2{Bo)  . 

Therefore 

(2.8)  \\u  —  Qu\\L2(B}  <  IMIl2(b)  +  II<5«IIl2(b)  —  ^  IMIl2(b)  —  CH  \u\Hi^  ■ 

Since  |'Ih|#i  <  C^H-  we  have 

\Qu\2hi{b)  <  max  {card  (A fi)}  E ies  (u) 


(2.9) 


<  C-H  e 


2  1  2 
z^£b  u  L2(n4)  -  Cjjg  u  L2(Bo) 


-  C BS H2  \u\Hi(Boj 


niL  ii  A 

°  S  a  H1  (Bq) 


The  estimates  above  are  for  the  set  B  around  the  edge  of  the  domain.  We  complete  the  proof 
by  making  similar  estimates  in  Qmt. 

Let  uB  be  an  extension  of  u  such  that  uB  =  u  on  H  and 


\ue\ Hl(nd)  —  C  (^)  lwlzfi(o)  ■ 

Let  B\  =  UjeAfi  and  Id  2i,  be  the  ball  around  B[.  Then  by  Assumption  1(c),  we  have  that 
diam  ( Bi )  <  CH. 


For  j  =  1, . . . ,  J,  define  Cj  =  JB.  uB  dx,  Uj  =  ue  —  cj.  Assumption  1.2  and  the  definition 
of  Q  imply  that  for  x  €  H  fl  Qint, 


Qu  (x)  =  Qui  (x)  +  Q. 


DOMAIN  DECOMPOSITION  PRECONDITIONER  FOR  GROUNDWATER  FLOW 


7 


Therefore 


U  ^  ^  '  ||ri 


We  use 


lie« 112 


i= 1 

J 


L2(ninnint) 


X  || Ui  +  Ci-  Q(ui  + 


i= 1 

J 


=  X]  llM*  “  QU 


1 1 2 


i=l 

J 


hlL2(Q.nf2mt) 


since  Qci  —  Ci 


<x 


|u"2 


i— 1 


+  ||Qu 


1 1 2 


*llL2(ninoint)  ■ 


< 


'*llL2(ninoint)  — 


EjeM  a,  (Ui)  ^|L(o)  -  Wj  («0 1  II^IIl2(o,))^ 


<  card  (TV*)  EjeM  («i)  ll^jll^n,)  <  CHd  ^jeAr.  a]  (Ui) 


<  CHdH  d  EjeM 


|«. 112 


<  C'Hm,"2 


*11  L2(f2j)  —  ^  II  u'*IIl2(_bj;)  j 


and  the  Poincare  inequality  to  obtain 

||W  —  Qu\\L2(pint}  <  2  Ej=l  ||«i||L2(B4)  <  CH2  Ej=l 


(*<) 


=  ^2  Ef=r  \uE\2m{Bi)  <  CH 2  zLi  M m(Bi)  <  CH2  \u\2m 


(Bi) 


l(Q)  • 


This  and  (2.8)  imply  (2.6). 

We  complete  the  proof  of  (2.5)  by  combining  (2.9)  with  the  estimate 

\Qu\iil(siint)  —  Ef=i  =  E*J=i  |Q  (“*  +  c*)liri(njnnin<) 


—  Ej=i  Qui  +  ci\' h1  (ninnint)  Ej=i  Qu 


<  Eli  EjeM  («i)  xI/j  Rl(n)  <  E/=i  card  (A/-,)  EjeJVi  («i)  l^l#i 


(«) 


—  ^  Hd  s  E*=i  EjeAi 


Kh"2 


—  C  hs  Ej= i  II^IIl2^) 


Il2(Oj)  -  ^  hs 


—  C hs^2  E?=i  "i  ~  C  d  Ej=i  —Cs  M#i(o) 


l„,l2 


This  completes  the  proof.  □ 
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LEMMA  2.3.  Under  Assumptions  1.1  and  1 .2,  for  every  finite  element  function  u  G  Vh,  there 
exists  a  decomposition  ut  G  Vu  such  that 

j 

(2.10)  u  =  Ylui 

i= 0 

(2-11)  X  lM*lfl'1(n)  —  C  +  ~f)  \u\ ■ 

i=0  \  0  / 

Proof  Define  Ih  to  be  the  fine  mesh  operator  Ih  '■  C  — >  Vh  by 

n 

h  («)  =  X«(*0ifc 

i=l 

where  {^}”=i  is  the  finite  element  basis  on  the  fine  mesh,  and  {xl}Tf1  are  the  fine  mesh  nodal 
points.  Let  u  be  partitioned  such  that 


u0  =  Qu  and  ut  =  Ih  (Ofiu  -  Qu )) , 

where  {$,}  is  a  partition  of  unity  such  that  9 fix)  =  1  for  x  G  9t  =  0  for  x  Qt.  and 

|V0i|  <  J-  Clearly,  by  construction, 

j 

X  “i  =  u- 

i= 0 

The  standard  arguments  [18]  imply  that 

(2.12)  X  lM*lfl'1(n)  —  ^  —  Qu\\ L2(n)  +  \u  ~  Qu\h1{q)  +  \Qu\H1(n)) 

i= o  vo  7 

We  complete  the  proof  by  applying  Lemma  2.2  to  estimate  \Qu\2Hl^  and  ||u  —  Qu\\'2L2^ny  □ 

If  we  let 

K°  =  C(1  +  j)  . 

we  complete  the  proof  of  Theorem  1.1. 

3.  Numerical  Results. 

3.1.  Laplace  Equation.  In  this  section  we  consider  the  simple  test  problem 
(3.1)  V2u.  =  0 

on  the  unit  square  [0, 1]  x  [0, 1]  with  zero  Dirichlet  boundary  conditions.  In  both  §  3.1.1  and  3.1.2 
we  use  the  function  identically  equal  to  one  on  the  mesh  as  the  initial  iterate  for  a  preconditioned 
conjugate  gradient  iteration.  We  terminated  the  iterations  when  the  residual  had  been  reduced  by  a 
factor  of  10“4. 

We  obtained  similar  results  using  GMRES  [17],  Bi-CGSTAB  [20],  and  TFQMR  [7]  in  the 
tests. 
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3.1.1.  Finite  Difference  Discretizations  with  Overlap.  The  experiments  in  this  section  were 
designed  to  measure  the  effects  of  overlap.  We  let  h  =  2~m  be  the  scale  of  the  fine  mesh  and  let 
the  overlap  o  be  the  is  the  nearest  integer  larger  than 

(3.2)  2moi  +  o0 

where  0\  is  the  overlap  that  depends  on  the  physical  domain,  and  oo  is  the  overlap  that  is  a  constant 
number  of  gridlines.  The  grid  is  an  n  x  n  mesh  where  n  =  2m  +  o.  In  this  way  we  can  subdivide 
the  region  into  4P  subdomains,  each  of  size  m  x  m,  where 

m  =  2m~p  +  0—1. 


The  scale  H  of  the  subdomains  is  2~p. 

We  considered  overlaps  of  1  (o0  =  1,  oi  =  0)  and  1%  (o0  =  0,  oi  =  .01).  We  discretized  with 
five  point  differences  and  use  the  function  identically  one  on  the  mesh  as  the  initial  iterate.  We 
used  the  MATLAB  software  associated  with  [12]  for  these  experiments  and  used  bilinear  coarse 
mesh  basis  functions.  These  do  not  form  a  partition  of  unity  as  a  piecewise  linear  set  would  and 
were  used  to  illustrate  the  flexibility  of  the  preconditioner. 

In  Table  3.1  we  tabulate  the  number  of  conjugate  gradient  iterations  needed  for  termination 
for  several  values  of  h  and  H,  with  an  overlap  of  1  and  1%.  As  the  theory  predicts,  the  iteration 
count  is  constant  as  h  and  H  are  reduced  simultaneously  for  the  overlap  of  1  and  declines  as  H  is 
reduced  with  the  overlap  of  1%. 


Table  3.1 

Iteration  Statistics  for  Two-level  Additive  Schwarz  with  CG 


Overlap  =  1 

Overlap  =1% 

H\h 

1/32  1/64  1/128  1/256 

1/32  1/64  1/128  1/256 

1/4 

12  16  21  28 

12  16  17  20 

1/8 

12  15  21  27 

12  15  17  19 

1/16 

8  12  16  21 

12  13  15 

1/32 

8  12  15 

10  11 

1/64 

8  12 

8 

1/128 

8 

10 


JENKINS,  KELLEY,  MILLER,  AND  KEES 


3.1.2.  Finite  Element  Discretization  with  Minimal  Overlap.  In  Table  3.2  we  report  similar 
results  for  a  piecewise  linear  finite  element  discretization  of  (3.1).  The  overlap  in  this  computation 
was  minimal. 


Table  3.2 

Finite  Element  Discretization 


H\h 

1/64 

1/128 

1/256 

1/4 

37 

51 

68 

1/8 

32 

44 

61 

1/16 

26 

36 

49 

1/32 

26 

37 

3.2.  Richards’  Equation  Results.  In  this  section  we  consider  a  test  problem  that  models 
groundwater  flow  in  a  partially  saturated,  homogeneous  soil.  The  model  equation  is  the  “head- 
based”  form  of  Richards’  equation, 


(3.3) 


06  Ss  dip 

dip  9S  dt, 


V  •  [KskrV  {ip  +  z)\ , 


where  ip  is  the  pressure  head,  9  is  volume  fraction  of  the  wetting  phase,  and  kT  is  the  relative 
permeability  of  the  wetting  phase.  The  remaining  terms  are  scalar  coefficients  given  in  Table  3.3, 
along  with  their  values  for  the  test  problem.  9  and  kr  are  functions  of  ip  given  by, 


(3.4) 

(3.5) 

(3.6) 


9={9S-  9r)  (l  +  |c#|n)  m  +  9r 


kr  —  (l  + 


—m/2 


j  aiPl"-1  (l  +  \a^P\n) 


ip  —  min  (ip,  0), 


where  m  —  1  —  1/n.  These  functions  are  derived  from  the  van  Genuchten  [21]  and  Mualem  [16] 
empirical  relations  among  pressure,  saturation,  and  relative  permeability. 

The  test  domain  is  the  unit  square  [0,  lm]  x  [0,  lm]  with  boundary  and  initial  conditions, 


ip(x,  0)  = 

0.0, 

Z  e  [0,  l] 

t  >  0 

ip(x,  1)  = 

0.1, 

X  e  [1/3, 2/3] 

t  >  0 

fOM)  = 

-1.0, 

x  G  [0,1/3)  U  (2/3,1] 

t  >  0 

= 

0.0, 

x  —  0, 1  z  G  [0, 1] 

t  >  0 

ip(x,z)  = 

#3 

x,  z  e  [0, 1]  x  [0, 1] 

t  =  0 

3.2.1.  Finite  Difference  Discretization  with  Minimal  Overlap.  We  discretized  equation  3.3 
by  applying  cell-centered  finite  differences  to  the  spatial  operator,  thereby  yielding  the  system  of 
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Table  3.3 

Richards  ’  equation  parameters 


Description 

Symbol 

Value 

Saturated  volume  fraction 

Os 

3.01  x  10”1 

Residual  volume  fraction 

0r 

9.30  x  10“2 

Specific  storage 

Ss 

1.00  x  10“6  (1/m) 

Hydraulic  conductivity 

I<s 

5.04  x  10°  (m/day) 

Mean  pore  size 

a 

5.47  x  10°  (1/m) 

Pore  size  uniformity 

n 

4.26  x  10° 

differential-algebraic  equations , 


(3.8) 


FijM,  f)  = 


d0_ 


+t(, 


hj 


dtpi. 


dt 


—z?  ~ 

~~KxI  “ 


^,3)  - 


1)] 


where  i  =  0, . . . ,  N  —  1,  j  =  0, . . . ,  N  —  1,  Az  =  Ax  =  1  /TV,  and 

(3.9)  ^!±jj  =  /2 

(3.10)  =  [(A".sA;r),J±i  +  /2 

The  semidiscrete  system  was  integrated  in  time  over  [0,  0.0149  days]  by  applying  the  fixed 
leading  coefficient  backward  difference  formulas  of  orders  one  to  five  [2, 11].  Order  and  step-size 
were  selected  via  local  truncation  error  estimates,  and  the  local  truncation  error  tolerance  was  set 
to  10/Ax2. 

At  a  given  step,  tn+ 1,  the  application  of  the  integration  method  yielded  a  nonlinear  system  of 
the  form, 

F[tn+i,ipn+i,  g(ipn+i)]  =  G(ifjn+ 1)  =  0 

where  g(y)  is  a  the  backward  difference  formula  for  dtp /dt.  We  solved  the  nonlinear  system  with 
an  inexact  Newton  iteration  that  terminated  when  the  2-norm  of  the  nonlinear  residual  was  reduced 
by  a  factor  of  10-5. 

At  each  Newton  iteration  we  obtained  the  Newton  step,  Sm+ 1 .  by  solving  the  linear  system, 


dG_ 

dtp 


(«i) 


■5ro+1  =  -G(*!n). 


with  scaled,  preconditioned,  BiCGstab.  The  scaling  was  obtained  from  the  integration  method’s 
weighted  root  mean  squared  norm.  Such  a  scaling  would,  in  real  applications,  allow  termination  of 
the  linear  iteration  according  to  tolerances  specified  by  the  integration  scheme  in  real  applications; 
however,  for  this  test  we  iterated  until  the  2-norm  of  the  true  linear  residual  was  reduced  by  a 
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Table  3.4 

Richards  ’  Equation  Iteration  Statistics 


factor  of  10  7  to  insure  that  errors  in  the  Newton  step  were  insignificant  with  respect  to  the  Newton 
iteration  and  integration. 

The  preconditioner  was  two-level  additive  Schwarz  with  the  coarse  grid  correction  determined 
from  the  restriction  and  interpolation  operators  in  §  1.3.  The  subdomains  had  the  minimal  overlap 
of  Ax  =  1  /N.  Table  3.4  gives  the  average  BiCGstab  iterations  per  Newton  iteration  for  two-level 
additive  Schwarz.  The  iteration  count  is  constant  as  H  and  h  are  reduced  simultaneously,  which  is 
consistent  with  the  predictions  of  the  theory. 
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